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The fractal properties of the total potential energy V as a function of time t are studied for a 
number of systems, including realistic models of proteins (PPT, BPTI and myoglobin). The frac- 
r— j tal dimension of V(i), characterized by the exponent 7, is almost independent of temperature and 

, increases with time, more slowly the larger the protein. Perhaps the most striking observation of 

this study is the apparent universality of the fractal dimension, which depends only weakly on the 
type of molecular system. We explain this behavior by assuming that fractality is caused by a 
self-generated dynamical noise, a consequence of intermode coupling due to anharmonicity. Global 
topological features of the potential energy landscape are found to have little effect on the observed 
fractal behavior. 
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I. INTRODUCTION 



The energy landscape of proteins near the folding temperature is often described as complexE'i That is, the 
landscape is expected to be rugged and composed of many minima and maxima, separated by barriers of varying 
heights, the result of frustration due to conflicting interactions between different segments af the protein. The complex 
nature of the underlying energy landscape has been demonstrated using ergodic measures^ and by exact enumeration 
^ of minima and transition states of a tetrapeptide.Q The goal of the present work is to investigate this energy landscape 
under a variety of conditions using tools from fractal analysis. The trajectories needed for analysis were generated 
' using "MOIL" , a molecular dynamics (MD) program which allows realistic simulations of small proteins using semi- 
empirical potentials!] Using MOIL, the total potential energy V(t) was calculated as a function of time, for several 
<^> proteins of choice (PPT - Pancreatic Polypeptide, BPTI - Bovine Pancreatic Trypsine Inhibitor and Myoglobin) at 
OC constant temperature. For comparison, we also simulated poly-alanine and NaCl. The V(t) curves so obtained under 

ON 

• a variety of choices of parameters, were then subjected to a fractal analysis. The main purpose was to correlate 
the resultant fractal dimensions with the topology of the energy landscape, and to understand the connection to the 
relevant physical parameters. Both MOIL and the method used to calculate the fractal dimension of the V(t) curves 
are described in Sec.||. 

There are several reasons for using a fractal approach in the analysis of the V(t) curves. Firstly, there is an intimate 
(— I connection between "roughness" and fractality: much is known about the -generation of rough surfaces using fractal 
Q , algorithms and conversely, on the analysis of rough surfaces in fractal terms.™ In this sense the notion of the roughness 
O ■ of the energy landscape of proteins naturally calls for a fractal analysis, which can quantify the degree of roughness. 
Secondly, fractal dimensions often show up as "universal exponents" , and are related to the critical exponents of phase 
transition theory.113 Conversely, there is a certain universality in protein dynamics, in the sense that most proteins 
exhibit a folding transition and multi-exponential relaxation times. This suggests a common feature in their energy 
landscapes. While the rapid transition to folded states in smaJLjproteins is best rationalized by the existence of a 
dominant native basin of attraction (NBA) or a folding funnelp'EJ the relaxation phenomenon might well be related 
to the roughness of the energy landscape, and hence to its fractal dimension. A "universal" fractal dimension would 
be a strong characterization of what it is that is common to different energy landscapes. Thirdly, establishing that 
the energy landscape of a protein is indeed fractal would allow one to relate the-subject of protein dynamics to a 
whole body of work on dynamics on fractals, e.g., diffusion on a fractal substrate. E3 
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The main problem in trying to characterize the roughness of a high-dimensional energy landscape, such as that 
exhibited by a protein, is the proper choice of a reaction coordinate. In this work we have chosen time for this purpose. 
The advantage is that time is a universal coordinate (unlike any of the spatial coordinates), and one can think of the 
protein potential energy dynamics in terms of this coordinate as a random walk taking place on a,rough substrate. 
Indeed, the rate of sampling of conformation space in complex many-body systems has been shownta to be a diffusive 
process in the energy space. The universality of the time-coordinate, however, can also be a drawback, in that it is 
not very sensitive to the global top olog y of the energy landscape, but rather to the local structure of minima. Indeed, 



as will be presented in detail in Sec. Ill, our results appear to contain information mostly on the role of anharmonicity 
in determining the fractal features of the V(t) curves. As such, the results reported here reflect more on the general 
properties of dynamics in anharmonic potentials than on protein dynamics per se. These and other issues will be 
discussed and analyzed in detail in Sec. IV. Conclusions and a summary are presented in SecjV]. Finally, to establish 



a simple reference model, an analytica 
presented in the appendix 



solution for the fractal dimension of the V(t) curve in the Rouse model is 



II. METHODS 
A. MOIL Simulations 

The MOIL program has been described extensively in the literature.! In a nutshell, it takes as input the x-ray 
structure of a given protein from the Brookhaven Protein Data Base (PDB), minimizes this structure, and runs 
molecular dynamics. MOIL makes it possible to simulate the energetics, dynamics and thermodynamics of systems 
that vary in size, up to several thousands of atoms. It accounts for a number of different energy contributions: van 
der Waals and electrostatic (non-covalent); bonds, angles, torsions and improper torsions (covalent): 

^electrostatic + ^ ^bond + ^ tangle + ^torsion + ^ Vim proper ■ (i) 

The sums are over all interacting pairs of particles, and the components are given as follows: 

Wvdw = -ir - -jsr- \ 2 ) 

a$ and hi are constants that depend on the type of atom i\ r.y is the distance between atoms i and j. 

^electrostatic — — ~, (3) 

erij 

where % is the charge on atom i and e is the dielectric constant. 

Vbond = kij (nj - (nj)) , (4) 

where kij is the bond force constant, (x) denotes the equilibrium value of coordinate x. 

Vanglc = Kijk (%fe — (^ijk)) 2 , (5) 

where Kijk is the angle force constant and 9^ is the angle between two sequential bonds, i.e., bonds connecting atoms 
i,3 and j, k. 

^torsion = a m ^klCOS m (<pijkl)- (6) 

m— 1 

Here k, I are four sequentially connected atoms; 4>ijki is the angle between the planes defined by atoms k and 
j, k, I; a m ^jki is a constant which depends on the types of atoms i, j, k, I and on the periodicity of the rotation m (see 
Ref.[| for details). Finally, 

^improper = A%'fcZ i^ijkl — (®ijkl)) 2 ■ (7) 

Here j, k, I are all connected to atom i and <&ijki is the angle between the planes defined by atoms k and j, k, I. 

In the present work, we considered only the time dependence of the total potential energy V due to all these terms. 
The time series for a number of systems were generated using numerical simulations employing cither the canonical 
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or microcanonical ensembles. As is customary, all simulations were preceded by a short initial annealing (heating) 
period. In the "equilibration mode" (used mainly for testing), the kinetic energies of all particles were rescaled in 
order to preserve the constant temperature (canonical ensemble). In the more important "no-equilibration mode" 
(default unless mentioned otherwise), no such rescaling was attempted and the system was allowed to reach thermal 
equilibrium at the specified temperature (microcanonical ensemble). 

The exponents characterizing the fractal nature of V(t) reported below derive from averages over groups of 25 runs 
differering only in the initial random velocity distribution. A number of tests were performed to ensure that our 
results are not artifacts of the particular integration schemes used. We checked that using more than 25 runs did 
not significantly affect our results statistically. Between groups of 25 but for fixed molecule-type, the runs differed, 
e.g., in temperature, addition of a water solvation shell, treatment of the terminus hydrogens, bond-freezing, and 
deletions of residues. Most runs started at a temperature of =1K and heated in IK per time-step increments to 
the final temperature Tf, at which the rest of the simulation was run. The total duration of the simulations, once 
Tj was reached, was usually 10 or 25psec, corresponding to 20000 or 50000 time-steps respectively, of 0.5fsec each. 
Exceptions to these conditions include some shorter and longer runs, with some variation in time-step size. They will 
be discussed in detail below. 



B. Fractal Analysis 

The V(t) curves were subject to a standard fractal analysis, with the purpose of calculating their fractal dimension 

7 = 2- a, (8) 
where a is the Holder exponent, defined in the next section. 



1. General Theory 

Given any set of points E = {(x,y)}, one can estimate their fractal diawpsion by a number of (in principle) 
equivalent procedures. A popular one is the so-called "Minkowski Sausage" ,ll3li-3 where one draws a circle of radius 
R centered at each point in E, and calculates the area A(R) of the union ("sausage") of all such circles. The fractal 
dimension is then found as 



ME) = lim 



2 Iog[A(r)] 



log(i?) 



slo P e{log(l/i?),log[A(i?)/i? 2 ]} (9) 



where the second (approximate) equality follows by assuming that the first equality holds for all R and multiplying it 
throughout by log(l/i?). This assumption is equivalent to a linear relationship between log(l/i?) and \og[A(R) / R 2 ], 
so that 7 is found as the slope of the corresponding log-log plot. If such a linear relationship is not found, it is taken 
in practice as an indication that E is not a fractal set. Conversely, i£_a.t least one order of magnitude of linearity is 
observed, the convention is to consider E a fractal set in that range.uJ It should be noted that rigorously speaking 
the curve of a function is generally not a self-similar fractal, but rather a self-affine one. The latter refers to those 
cases where the axes have different units, so that the scaling might be different in the two directions. Indeed, a 
(deterministic) self-affine curve satisfies the following scaling relation: 

h(t) = b- a h(bt), (10) 

which holds only on average for random processes. Xku "Holder" , or self-affinity exponent a, is 2—7 [Eq.(||)] ifXhe curve 
is subjected to the Minkowski analysis of Eq.(|9|).ESEZl An example is the Weierstrass-Mandelbrot function:li3 h(t) = 
b~ na (l — cosb n t), which satisfies Eq. (|l0|) and whose fractal dimension is 2-<x An equivalent characterization 
of statistically self-affine fractals is through the correlation function, which satisfies: 

C h (t) = ([h(t + t )-h(t a )] 2 )~t 2a . (11) 

This relation is expected to hold for durations shorter than a typical correlation time t Cl found from: 

fdttT(t) , x 

c ~ JdtT(t) ' 1 ' 



where T(t) = (h(t )h(t + t )) - (h(t ) 2 ). 
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2. Analysis of Time Series 

In the present case the (self-affinc) set under consideration is: E = {(t, V(t))} t , and the fractal analysis we applied is 
a numerically well-behaved version of Eq. (0): the "e- variation method" . A full description of the method is beyond the 
scope of this paper; for full details see Ref.|18|. Briefly Dubuc et al£3 introduced a method particularly well suited for 
the evaluation of the self- affinity exponent. They demonstrated that their method has the most stable local exponent in 
comparison to a variety of other methods, such as box-counting and power-spectrum. In their notation, the e-variation 
of a function / is V(e, /) = J Q v(x, e)dx, where the e oscillation v(x, e) is: v(x, e) = sup x / eR t x -\ f[%') — inf x 'eR c (x) f( x ')i 
and where R e (x) — {s G [0,1] : |sc — s| < e}. The corresponding log-log plot for the calculation of the self-affinity 
exponent is [log(l/e), log(V(e, /)/e 2 )] , with the exponent given by 2-slope. We applied this analysis in order to 
extract the self-affinity exponent (which we refer to simply as the "fractal dimension" and denote by 7) of our V(t) 
curves. 



III. RESULTS 

In this section we will present the findings from our numerical simulations. A discussion will be postponed to the 
next section. We will start by presenting a sample of "raw" simulation results, then move on to a detailed presentation 
of each of the three proteins, compare them amongst each other, perform a comparison with several test cases, estimate 
the effect of time-scales of different processes, and finally study the role of anharmonicity vs harmonicity. All our 
MD simulations started out from the folded state (as determined by X-ray), which we subsequently minimized using 
the Powell algorithmic The protein was then heated from IK in increments of IK at every time-step, with a random 
initial velocity distribution, until the desired final temperature was reached. Fractal analysis was performed on data 
from the constant temperature period only. Dynamics on all three proteins was run for 25psec, with a time-step of 
0.5fsec (with a few exceptions, to be detailed below). Fractal analysis was performed on the full 25psec trajectories 
and on the first lOpsec of each such trajectory. This process was repeated for a wide selection of temperatures in the 
range 1K-800K, with a focus on the room temperature regime. 

A. "Raw" Results 

A graph of a typical V(t) curve is displayed in Fig.[l], for BPTI at 310K. We only present the time series after the 
temperature is stabilized, i.e., the first 0.4psec during sthich annealing from IK occurs are omitted. The noisy data 
are strikingly similar to "fractional Brownian motion" ,113 and clearly suggest the use of fractal- type analysis. This is 
done in FigJ|, where the corresponding e-variation analysis is shown. A straight line can be fitted with confidence 
to the e-variation data over some 1.5 orders of magnitude (decades). The linear regression coefficient in this case is: 
1Z 2 = 0.998, where 1 indicates a perfect straight line. The fractal dimension (slope in Fig.^J) is: 7 = 1.762 ± 0.004, 
where the uncertainty is due to the linear regression. All log-log plots look strikingly similar to the one in Fig.||, and 
have a similar range of linearity. The 1Z 2 coefficient rarely went below 0.99. Fig.|| presents all 25 fractal dimensions 
7 obtained from the runs for BPTI under the abovementioned conditions. The dashed line is the average fractal 
dimension, which in this case yields: 7= 1.75 ± 0.02. The uncertainty in this number reflects both the individual 
linear regression errors on each data point, as well the standard deviation about the average. The procedure illustrated 
in Figs.|l[j^ was repeated throughout this work. Thus every 7 data point in the figures to follow is the result of an 
average over 25 fractal dimensions obtained from independent runs differing only in the initial velocity distribution. 

B. Myoglobin 

Fig.|] shows the results of our simulations of myoglobin. Every 7 point was obtained as discussed above, so the error 
bars are due to both the averaging over the 25 data points and the linear regression error on each of these points. 
The dashed (solid) line connects the 7 values resulting from the 25psec (lOpsec) runs. The inset shows the room 
temperature results. In the lOpsec case 7 appears to have a shallow minimum around 10K, which is however absent in 
the full 25psec simulations. The error bars are in fact too large to confidently discuss a trend in the curves. However, 
it is clear that the fractal dimension increases with the simulation (or observation) time t: on average 7~ 1.67 for 
lOpsec and rises to 7« 1.72 for 25psec. This behavior is seen in all our results ahead (see Figs. || and g). 

As is clear from the inset, the fractal exponent 7 is essentially constant around room temperature, apart from the 
rather sharp drop at 330K. 
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The robustness of the results was tested by running the simulations at 300K and 400K with the heme group deleted. 
The removal of the heme group makes myoglobin considerably less compact, causing it to unfold rather easily. Hence 
we expected a significant change in the underlying dynamics. However, Fig. [I] shows (squares and circles) that the 
influence of heme deletion is very minor, and results in a small increase in 7. 



C. BPTI 

The results of our simulations of BPTI are shown in Fig.|5|. The overall behavior is very similar to that of myoglobin. 
The differences are (1) overall the fractal dimensions are higher; (2) the minimum is more pronounced (beyond 
statistical error), visible also in the 25psec data, and has shifted to 50K. 

The increase in 7 with increasing observation time for the longer simulation is unmistakable for BPTI as well. 

We also performed a lOpsec simulation including a water solvation shell (diamond). Clearly, this has no significant 
effect on the results. 

As a consistency test, we calculated 7 for the complementary 10-25psec of the simulation (triangle). The result is 
a 7 value in between the lOpsec and 25psec values. 

Finally, the dip observed for myoglobin can be seen here as well as in the 25psec simulation, at T=320K. 



D. PPT 



Our PPT simulation results appear in Fig]§ The general qualitative similarity to BPTI and myoglobin is noticeable, 
namely the increase in overall fractal dimension with simulation time, and the shallow minimum at around 50K. Several 
additional details are noteworthy: (1) We included a 50psec long simulation at IK (same 0.5fsec time st ep), and the 
increase in 7 persists (full circle). (2) We also performed simulations in the "equilibration mode" (see Sec. II A). These 
are indicated as squares and can be seen most clearly in the inset. The results are virtually identical to the usual 
"no-equilibration mode" , demonstrating the agreement between the canonical and microcanonical ensembles for large 
(> 100 atoms) systems. The only exception is the point at 800K, where the "no-equilibration" value of 7 drops rather 
sharply. (3) A water solvation shell (diamonds, also in the equilibration mode) again does not produce statistically 
significant results. 



E. Comparison Among Myoglobin, BPTI and PPT 



Figs.0 and I compare the results for the three proteins in the 10 and 25psec cases, respectively. These figures merely 
recapitulate the results from FigsJ^-^ and show no new data. The main points to notice are the striking similarity 
between PPT and BPTI, compared to the overall lower fractal dimension values for myoglobin. This holds in both 
the 10 and 25psec cases. 



F. PPT Under a Variety of Conditions 



We performed further tests of the robustness of our results by running lOpsec simulations of PPT at 300K, under 
different conditions, as shown in Fig.^|. There are 11 data points in this figure, and the 7 values are arranged in 
increasing order. Let us briefly describe the results, noting that case 7 is the "standard" one (i.e., the result shown in 
Fig|j). 

The case "no a- helix" refers to a simulation in which the entire PPT a- helix was deleted, i.e., only the GLY PRO 
SER GLN PRO THR TYR PRO GLY ASP ASP groups remain. Under these conditions a significantly smaller 7 value 
is obtained. On the other hand data point 5 is the opposite case, where only the a-helix remains. The corresponding 
7 is much closer to the reference value of data point 7. This shows that the a-helix is by and large responsible for the 
observed fractal dimension. 

Data point 2 corresponds to a freezing of all bond vibrations. The remaining motions are therefore proper and 
improper rotations and center of mass translations. Since the deviation from the reference 7 value is small, it appears 
that the rotations are largely responsible for the observed fractal dimension. 

Data point 3 is the result when a different integrator is used, which appears to have only a minor effect. This test 
was performed to check whether the results are a consequence of numerical instabilities. 
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Data point 4 is another deletion experiment, this time of just the SER residue. This has a very small effect which, 
curiously enough, is slightly larger than leaving only the a-helix, as indicated in case 5 (although due to the large 
error bars in fact this effect cannot be taken too seriously). 

To calculate the 7 corresponding to data point 6, the terminal hydrogens were treated as "extended": unlike in all 
other cases, the hydrogens were not treated as atoms but were "absorbed" into the carbons. This modifies the carbon 
mass and radius to appropriate effective values. The effect, however is negligible. 

Data points 8 and 10 relate to using the equilibration mode and adding a water solvation shell respectively, which 
(as shown also in Fig.|J) do not have a significant effect either. 

Data point 9 reports the result from a single 250psec long trajectory: the trajectory was divided up into 10 
segments of 25psec each, and the 7's were averaged. Data point 11 is due to a lOOpsec long trajectory which, however, 
employed a time step of 5fsec (10 times larger than the usual time step). In accordance with the observation of a 
fractal dimension which increases over time, the resulting average 7 of data point 9, as well as that of data point 11, 
are larger than the reference value. 

G. PPT vs Test Systems: Poly-alanine and NaCl 

In order to estimate to what extent our results reflects properties unique to realistic proteins, we present calculations 
for two test systems: poly (16)-alanine and a 40-atom cluster of NaCl. Poly-alanine folds into an a-helix, so it has the 
crucial part of the protein topology (recall points 1 and 5 in Fig.^|), but the underlying energy landscape is not expected 
to be rugged. This is becuase all the residues are identical so there is no frustration due to conflicting interactions 
between different segments of the polymer. (Nevertheless each residue has internal conflicting interactions, giving rise 
to some ruggedness). The ground state of the non-polymeric NaCl cluster is easily accessible!! i.e., starting from a 
disordered state, it finds it crystalline ground state well before the end of our lOpsec simulation. Thus it has neither 
the protein topology, nor its rugged energy landscape. The results of our simulations are shown in Fig.|l^. The 10 and 
25psec runs of 16-alanine and NaCl at 250K, as well as lOpsec at 300K, and lOpsec for NaCl alone at IK. For clarity 
of presentation the error bars have been removed. However, it should be remarked that when the error bars are taken 
into account the results for a given temperature overlap. Indeed, it is striking that the results for poly-alanine and 
NaCl do not differ significantly from those for PPT. In fact for 250K, NaCl and PPT agree almost exactly, and the 
agreement for lOpsec at both 250K and 300K is quite close. Only at IK does there seem to be a significant difference 
between NaCl and PPT. 16-alanine has lower 7 values than PPT, but not very much so. Barring accidental agreement, 
these results clearly indicate that neither the protein topology, nor the rugged energy landscape are essential in order 
to obtain fractal V(t) curves, with fractal dimension close to that of realistic proteins. 

H. The Role of Observation Time: Short Simulations of PPT 

In order to systematically investigate the role of the total duration of the simulation and the correspondingly 
participating physical processes, we ran simulations of 20000 time steps of PPT at 300K, with progressively smaller 
time steps. This is different from the bulk of our simulations where a constant time step of 0.5fsec was employed, 
so care should be exercised in the comparison. The main point is that by employing a smaller time step, previously 
inaccessible fast motions are now observable. The results are shown in Fig.jll]. It is remarkable that the fractal 
dimension obeys an almost perfect logarithmic law over three orders of magnitudes: 

7(r) = 1.60 + 0.060 ln(r), (13) 

for O.lpsec < r < lOOpsec. The regression coefficient is_7^. 2 = 0.999. Fig.[n] thus demonstrates decisively that 7 is 
indeed strongly dependent on the total observation time.Q 

Below t = lOOfsec (to the left of the arrows in Fig.|ll|) the behavior of 7 appears rather erratic and we believe this 
is due to numerical error. Indeed, the fastest motions in the protein are the OH or CH stretches, with a period of 
about lOfsec. Thus r = lOfsec would most likely not probe even a single period and can confidently be regarded as 
the lower physical limit of our simulations. 

We turn next to an analysis and discussion of our results. 

IV. DISCUSSION 

The most striking observation presented in the previous section is the rise of the fractal dimension with time. It 
suggests (at least) two likely hypotheses for its explanation: As time evolves 
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1 . there is increased access to local minima and "fine-structure" of the potential energy landscape; 

2. low frequency vibrations and energy transfer between modes become increasingly activated. 

The first of these attributes the growth in 7 over time to the increased roughness of the sampled potential energy 
landscape. In other words, the protein "visits" more and more local minima (or basins of attraction) as its samples 
its conformational space. This process is consistent with an increase in fractal dimension since the fractal potential 
energy landscape is not sampled all at onceJpy the protein: instead the inner details are only progressively revealed. 
However, as shown by Czerminski and Elber,a activated transitions between similar conformational structures confined 
to a small set of minima (separated by barriers on the order of ksT) can occur at room temperature on a time scale of 
several picoseconds. A much larger time scale is required for such transitions at temperatures on the order of IK. The 
insensitivity of 7 to temperature therefore appears to rule out the first hypothesis: at the low temperature end one 
would expect the protein to be confined to a single minimum within the time scale of our simulations, and the fractal 
dimension to be smaller as a consequence. Thus while increased access to local minima is consistent with the growth 
of 7 with time, it is inconsistent with the lack of temperature dependence. We conclude that the fractal dimension 
measured in our simulations is not a good probe of the global potential energy landscape topology. 

This brings us to the second hypothesis which, as will be argued shortly, is consistent with both the observed time 
and temperature dependencies. According to this hypothesis, the primary factor determining the fractal dimension of 
the potential energy curves is the self-generated dynamical noise, namely the very broad frequency spectrum appearing 
in the dynamics. This spectrum and the resulting self-generated noise is a consequence of the anharmonicity in the 
potential, as is demonstrated convincingly in Fig.|l^. Shown in this figure, along with the full lOpsec simulation of 
PPT at 300K, are the results of the same simulations with all anharmonic interactions switched off: terms K,dWj 

^electrostatic 

and Vtorsion are missing from the total potential energy V of Eq.([IJ). The effect is dramatic: there is a 
strong decrease of 7 as the temperature is lowered. This shows that anharmonicity is of paramount importance in 
determining 7, at least up to 100K. 

The mechanism we propose is thus as follows: V(t) is a sum over many anharmonic terms, the number of which is 
proportional to the number of atoms in the protein. These anharmonic terms create overtones, which effectively fill the 
frequency spectrum and create a quasi-continuum. Now, on the one hand, as time evolves (at constant temperature) 
both very low and very high (previously inactive) frequency components become activated: the low frequencies are 
only sampled on long times, whereas the high frequencies correspond to highly energetic modes which must await 
energy transfer before equilibration sets in. On the other hand, as temperature increases (at constant r) intermode 
coupling causing energy transfer from rigid modes to soft modes is facilitated, and more frequency components are 
activated. Thus in both cases one expects an increase in dynamical noise, with a corresponding increase in the fractal 
dimension of the signal. This accounts for the behavior seen in Figs.^|-[| Unfortunately we have not been able to 
study whether the fractal dimension converges on a limiting value upon full equilibration (at least as a function of 
temperature this does not seem to be the case - see Figs.[j] and |^), since this requires prohibitively long simulations.!] 

The behavior of the harmonic case observed in Fig.|l2| might seem curious in light of these arguments: after all in 
the purely harmonic case there is no intermode coupling, nor any dynamical noise. However, we recall that due to 
the use of non-Cartesian coordinates [see Eqs.Q), (JbJ) and (Jt|)], a certain residual degree of anharmonicity remains 
with respect to the Cartesian components. At very low temperatures vibrations and rotations occur very close to 
their equilibrium values and the anharmonic effects are barely noticeable. Only as the temperature rises do these 
vibrations and rotations start to deviate significantly from their equilibrium values and gain anharmonicity. This 
explains the rise of 7 in the "harmonic" case. In contrast, the other terms [Eqs.(pl), (||) and (^j)] make a strong 
anharmonic contribution at all temperatures, which accounts for the relative constancy of 7 in this case. 

The different behavior of PPT and BPTI vs myoglobin can be attributed to the size of these proteins. Myoglobin has 
156 residues, whereas PPT and BPTI only have 41 and 60 residues respectively. Therefore at any given temperature 
it takes longer to activate the overtones that give rise to the self-generated noise in myoglobin than in either PPT or 
BPTI, which explains the lower fractal dimension of myoglobin in Figs|7]]|. The proximity in size between PPT and 
BPTI, on the other hand, is the reason for their quantitatively similar behavior. A related argument makes clear why 
deletion of the heavy heme group from myoglobin yields a higher 7 (Fig.^J): without heme the energy redistribution 
among the modes occurs faster. 

Also the curious similarity between NaCl and PPT can be attributed to anharmonic effects: NaCl, with its purely 
Coulombic interactions, is of course strongly anharmonic, even more so than PPT. On the other hand poly-alanine is 
more harmonic, so as expected it has lower 7 values. 

It is thus seen that the "anharmonicity" hypothesis is capable of accounting for all of the salient features present 
in our results for the fractal dimension of the V(t) curves. The role of multiple minima and ruggedness of the energy 
landscape seems to be more restricted: it is apparently not the global topology which is measured by 7. Nevertheless, 
it should be noted that different multiple minima will generally lead to a denser filling of the frequency spectrum or 
to the appearance of bands, which may overlap and thus induce additional intermode coupling. 
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An important remaining issue concerns the value of 7 observed in our simulations: definitely greater than 1.5 (except 
perhaps for the r = lOOfsec simulation, Fig.|Tl|). A priori it is not clear what the lower limit of 7 for polymer-like 
systems subject to the usual dynamics (like Newtonian or Langevin) should be. A plausible lower limit of 1.5 is 
derived in the appendix for the Rouse model. The Rouse model considers a purely harmonic set of beads under the 
influence of a Gaussian random force. In our simulations (in the "no-equilibration" - microcanonical mode) there 
is no randomness except for numerical truncation errors and "randomness" due to chaos caused by the anharmonic 
terms. However, we find 7~ 1.5 as the result of our "harmonic" simulations at very low temperature (Fig.|l2|). As 
argued above, this actually corresponds to the low-anharmonicity limit, so is in agreement with the result for the 
Rouse model, provided the truncation errors and/or the chaotic "randomness" can be modelled as a Gaussian random 
force acting on the protein (this is a speculation). Increasing the degree of anharmonicity can only increase the fractal 
dimension of the signal, by the mechanisms discussed before. These arguments then show that the Rouse model result 
could thus serve as a lower bound to the 7 values found in our simulations. 



V. SUMMARY AND CONCLUSIONS 



In this work we have studied the fractal properties of the total potential energy V as a function of time t for 
several realistic proteins. We performed simulations at a wide range of temperatures and a number of time scales. 
Our main conclusion is that the fractal behavior seems to be caused by the self-generated dynamical noise due to 
intermode coupling and equilibration promoted by anharmonic effects. There are universal aspects to this behavior, 
as exemplified by the similarity between totally different systems such as PPT and NaCl. Thus it cannot serve to 
characterize individual proteins and their energy landscape, but is rather a property common to all sufficiently large 
and anharmonic systems. We believe, however, that probing the fractal aspects of protein potential energy landscapes 
is a useful approach to the characterization of the ruggedness of these landscapes. To achieve this, it would be very 
interesting to examine V as a function of appropriate spatial coordinates, rather than t. The question remains whether 
there are other ways leading to useful geometrical characterizations of the energy landscape, e.g., identification of 
multiple minima and maxima by real and imaginary frequencies a 

Mpte added: While this paper was in the final stages of preparation we became aware of an article by Garcia et 
al£3 The authors concluded that the time dependence of the mean-square displacement of crambin around the crystal 
structure is described by an effective fractal exponent 7 ~ 1.60 (for t long). This result was interpreted as indicative 
of multi-basin dynamics which would be typical in a rugged energy landscape J3 
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APPENDIX 



AN ANALYTICALLY SOLVABLE EXAMPLE: THE ROUSE MODEL 



It is useful to have a simple reference model in the background to compare the numerical results of Sec. Ill with. 
To this end, we consider the Rouse Model (RM), which is standard in polymer theory. It describes a polymer as a 
set of N beads located at (Ri, R2, Rjv) connected along a chain, where each bead (n) undergoes Brownian motion 
under the influence of random force f„(i). It should be remarked that this is essentially different from our simulations, 
in which no randomness (except for numerical truncation errors and chaotic "randomness" due to anharmonicity) is 
present. Nevertheless, it useful to study the RM in order to see analytically how a fractal V{t) can come about. 

The excluded volume and-hydrodynamic interactions are disregarded in the RM. The dynamics is generally described 
by the Langevin equation:E3 

§-n n{ t) = £ H nm • ( U t) - ™) + I fcs r£ * • n nm . (Ai) 
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In the RM the mobility tensor is assumed to be isotropic: 



and the interaction potential is taken to be harmonic: 

k N 

^ = -^(R„-R„_ 1 ) 2 , (A3) 



2 

n=2 



so that Eg. ([Al|) reduces to 



C-^ = -fe(2R„-Rn+i-Rn-i)+fn (n = 2, 3, N - 1) (A4) 

C^ 1 = -fc(Ri - R 2 ) + fi (A5) 

c dRjv = _ fe(Rjv _ Rjv _ l)+fiV . (A6) 
at 

The force and the potential define a time scale through 

r-f <A7) 

It is convenient to consider the last set of equations in the limit where n is taken as a continuous variable: 

<9R„ <9 2 R ra 

C^ = fc ^TT+ f w (A8) 



dt dn 



0. 



dn dn 

The random forces are assumed to be Gaussian, i.e., delta-function uncorrelated: 



&.(*)> =0 (A9) 
<f„(t)f m (t')> = 2Cfc B T^(n - m)5(t - t')- 



In the continuum limit the potential energy becomes: 



2 J \ dn 
The last three sets of equations define the continuum RM 



1. Normal Coordinates 



U=U dn(^) . (AH)) 



We are interested in calculating the fractal scaling of the potential energy U (t) defined in Eq. ( AlOQ . For this purpose 



it is convenient to transform to normal coordinates. It can be shown that in terms of the coordinates 



i r 

X p =— / dn cos{pnn/N)R n (t) p = 0,1,2,... (All) 
Jo 



the Langevin equation (AE) becomes 



Cp ^ = -k P X p + f p (A12) 



where 
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and the forces satisfy 



Co = N( (A13) 

( P = 2N( for p=l,2,... (A14) 

2ir 2 kp 2 , s 

k p = ^ L for P = 0,1,2,... (A15) 



(f p (t)) = (A16) 



(f p (i)f g (0) = 2( p k B T5 P q5(t - t'). (A17) 
In this representation the random forces as well as the motions of the X p 's are independent, so that the motion of the 



polymer is decomposed into independent modes. Indeed, it is easily verified that the formal solution to Eq.(A12) is 

X p (t) = ^f t dt'e-C-*')/^*'), (A18) 



where 



r p = ^ , n = ^ = ^-r. (A19) 
p z ki tt z 



Finally, the inverse transform of Eq. ( All 



is 



R„ = X + 2 ^2 X p cos(pirn/N). (A20) 
P =i 



2. Scaling of the Potential Energy 

We will now present a calculation of the fractal scaling of the potential energy in the RM, as defined in Eq.(|ll|):0 



Cu(t) = ([U(t)-U(0)f). 
Inserting the expression for R n into Eq.( Aioj) and expanding yields: 



U = 2k (-) P9 x p(*) ' X <?W / dn sm{pTm/N)sm(qim/N). 

p,q=l •'° 



The integral evaluates to ^-S pq , so that: 



2 N 



P =i 



Thus the fluctuation can be written as: 



Cuit) 



kir 2 
kix 2 



N 



Yp 2 (X^(i)-X^(O)) 



.p=l 



^ ) £ P V [<X^(t)X»(t)) - 2<X*(i)X*(0)) + (Xj(0)X»(0))] 



p,q=\ 



Inserting the formal solution for the normal coordinates, Eq.( A18| ), the averages can be calculated from: 



(A21) 



(A22) 



(A23) 



(A24) 
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<xj(t)x»(f )> = ^ /_* 



dti / dt 2 I dts / eft 4 x 

Sp >/ — oo •/ — oc J — oc J —oo 
-(a*-*i-*»)/r, e -(«'-t8-t4)/r e ^ (tl j.^ (ta)f?(t3) .^ (t4) j 



(A25) 



To proceed the averages over the forces must be evaluated. This can be done with the help of Wick's theorem:E3 

(A26) 



{X ni X n2 ■ ■ ■ X n2p ) S ' (imi^n^) ' ' ' {x m2p l X m2p ') 

all pairings 



where the :r n 's are Gaussian random variables. Employing the orthogonality relations from Eq.(Al7) thus yields: 

(f P (ti) • Ufa) f,(t 3 ) • f q (U)) = 
(2k B T) 2 [( p ( q 6(h - t 2 )S(t 3 - U) + CpS 2 pq S(h - t 3 )6(t 2 - U) + Qpl^ih - U)5{t 2 - t 3 )] (A27) 

In addition we need the rule 

dh f(t - t 2 )S(t 2 - t x ) = f(t ~ h)G(t ~ h) (A28) 



whe re Q (x) is the Heavyside step- function. Combining these results, the delta functions simplify the expression in 
Eq.(|A2|) as follows: 



X£(i)X2(f )) = (2k B TY 



J_ I dtl e -(t-ti)(i/r p+1 /n) f dt3 e -(f-t3)(l/r,+l/n) + 



o^P9 



min(i,r, ) 



{2k B Tf 



T p T q \ 1 



ip T 'g / ^p^i? S; 



+ 9A ZP.„2(2min(t,t')-t-t')/T P 
P 



(A29) 



Inserting this into Eq.( A25 ) yields: 



2 2 



\ 2 1 T 2 



Tp + T 9 / CpC<? Cp 



4 I| (l _ e 2*/r P 



'A"/ 



(A30) 



Using the definitions of the various parameters appearing here the final result becomes: 

N , N 



Cu(t) = (k B Tf 



.p^q— 1 p— 1 



1 — e r ^ 



(A31) 



The time dependence of Crjit) for large but finite TV arises from the term Ylp=i ex P( — 2ir 2 tp 2 /tN 2 ). To derive an 
expression for 7, it is useful to approximate the above sums by integrals. Defining x = p/N, y = q/N , we have: 



C v (t) « (k B T) 2 



47V H / dx I dy 



N 2 Ji/N Ji/N x 2 + y 2 J 1/N 



4 / dxe~ 2 * i x 



(A32) 



The first integral is easily evaluated after a transformation to polar coordinates, and the second corresponds to the 
error-function (in the second integral the error in taking 1 /N = is negligible) : 



Cu(t) « (k B Tf 



' Ar „ IniV „ ^crfVe 
47V + 2.— -27^ 



, 9 = 2ttH/t. 



(A33) 



In accordance with the preceding comments we next take the t -C r limit. Then, using erf(a;) rs -^(x— + — •■■) 
to the first non-trivial (third) order, we can write: 
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4^ + 2^-4(1-^) 



(k B T) 



,8tt 2 t 



2(2- 7 ) 



(A34) 



where in the last line the time-independent terms were ignored, and ultimately 7 = 3/2. Thus the Rouse model has a 
potential energy fluctuations which have a fractal dimension of 3/2, just like a classical, fully correlated random walk. 
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FIG. 1. Total potential energy as a function of time for a single run of BPTI at 310K. Annealing from IK occurs during the 
first 0.4 psec, which are not shown. The highly irregular fluctuations are reminiscent of fractional Brownian noise. 




FIG. 2. e-variationliS result for the curve shown in Fig^]. A straight line can be observed for about 1.5 decades. The bending 
at high e values is due to finite size effects. 
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FIG. 3. All 25 fractal dimensions calculated for BPTI at 310K. Individual error bars are due to linear regression. The dashed 
line is the average fractal dimension. 
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FIG. 4. Fractal dimension as a function of temperature for myoglobin. Lines connecting data points are guides to the eye 
only. 




FIG. 5. Fractal dimension as a function of temperature for BPTI. 
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FIG. 6. Fractal dimension as a function of temperature for PPT. 
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FIG. 7. Comparison among myoglobin, BPTI and PPT for lOpsec simulations. 
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FIG. 8. Comparison among myoglobin, BPTI and PPT for 25psec simulations. 
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FIG. 9. A comparison of results for PPT at 300K, lOpsec, under a variety of conditions. 
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FIG. 10. A comparison between PPT, 16-alanine and NaCI (40 atoms). Error bars have been omitted for clarity and are 
replaced by small dots. 
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FIG. 11. 20000 time step simulations of PPT at 300K with varying total duration (x-axis). Results to the left of the arrows 
are in the unphysical range. 
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FIG. 12. lOpsec simulation results of PPT with the full potential (solid line) and just the harmonic part (dashed line). 



